cell_dict={}
median_dict={}
min_x = 0
max_x = 9
mix_number=2
max_mix=20
EB_cell_number=20000
total_cycles=10000
nFinalCells = 100000
minCellsPerCelltype = 100
RarestCelltypeRate = 0.09
#what is the minimum number of cells per genotype given RarestCelltypeRate and minCellsPerCelltype required?
# RarestCelltypeRate(x100) : 100 = minCellsPerCelltype : X
MinCellsPerGenotype = (100*minCellsPerCelltype)/(RarestCelltypeRate*100)
MinCellsPerGenotype_Rate = MinCellsPerGenotype/nFinalCells
print(MinCellsPerGenotype_Rate)
x, probabilities = empirical_dist(GrowthRate["D0_to_D4_growthRate"].tolist(), min_x, max_x)
x1, probabilities1 = empirical_dist(GrowthRate["D4_to_D10_growthRate"].tolist(), min_x, max_x)
k=0
while mix_number <= max_mix:
cell_dict[mix_number]={}
for i in range(1,total_cycles):
generated_data = np.random.choice(x, size=mix_number, p=probabilities)
generated_data_2 = np.random.choice(x1, size=mix_number, p=probabilities1)
if sum(generated_data)>mix_number:
k+=1
generated_cell=(EB_cell_number/mix_number)*generated_data
# if sum(generated_cell*generated_data_2)>sum(generated_cell):
# cell_dict[mix_number][k]=generated_cell*generated_data_2
cell_dict[mix_number][k]=generated_cell
simulated_exp=pd.DataFrame(cell_dict[mix_number])
plt.clf()
plot_data=(simulated_exp.div(simulated_exp.sum(axis=0), axis=1)>MinCellsPerGenotype_Rate).sum(axis=0)
testplot=sns.histplot(plot_data,discrete=True,stat='probability')
testplot.set(xlabel='Number of analysable genotypes', ylabel='Frequency')
testplot.set(xticks=range(0,15,1))
plt.title('Number of genotypes mixed: '+str(mix_number))
# plt.show()
median_dict[mix_number]=[plot_data.mean()-plot_data.std(),plot_data.mean(),plot_data.mean()+plot_data.std()]
mix_number+=1
print(plot_data.std())
plt.clf()
data_median=pd.DataFrame(median_dict).T
data_median['x']=data_median.index
#%%
# PLot empirical distribution
sns.set_style("whitegrid")
sns.set_style("ticks")
fig,ax = plt.subplots(ncols=1,nrows=1, figsize=(7, 5),dpi=100 )
sns.set_style("white")
#plt.plot(data_median["x"], data_median[1], '-', linewidth = 3, label = "Median")
sns.lineplot(data=data_median,x="x",y=1,color="#FDA63A",
linewidth=3, ax = ax)
sns.lineplot(data=data_median,x="x",y=0,color="blue",
linewidth=1.5, ax = ax)
sns.lineplot(data=data_median,x="x",y=2,color="blue",
linewidth=1.5, ax = ax)
plt.fill_between(data_median["x"], data_median[0], data_median[2] ,
alpha=0.2, color='blue',linewidth=0.0)
#plt.xlabel('pseudotime', size=30)
ax.set_ylim([0,
round(data_median[2].max()+1) ])
ax.set_xlim([0, data_median["x"].max()])
ax.set_xticks([0, 2,5,10,20])
ax.set_xlabel('Mixed cell lines') # Y label
ax.set_ylabel('Recovered cell lines') # Y label
sns.despine(fig=ax.figure)
fig.savefig(figDir+"/MinGenotypeRate0.05.Simulation.RarestCT:{}.svg".format(RarestCelltypeRate), format='svg', bbox_inches='tight')
